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Abstract 

Dielectric  Barrier  Discharge  (DBD)  type  devices,  when  used  as  plasma  actuators, 
have  shown  significant  promise  for  use  in  many  aeronautical  applications. 
Experimentally,  DBD  actuator  devices  have  been  shown  to  induce  motion  in  initially  still 
air,  and  to  cause  re-attachment  of  air  flow  over  a  wing  surface  at  a  high  angle  of  attack. 
This  thesis  explores  the  numerical  simulation  of  the  DBD  device  in  both  a  ID  and  2D 
environment.  Using  well  established  fluid  equation  techniques,  along  with  the 
appropriate  approximations  for  the  regime  under  which  these  devices  will  be  operating, 
computational  results  for  various  conditions  and  geometries  are  explored.  In  order  to 
validate  the  code,  results  are  compared  to  analytic  or  experimental  data  whenever 
possible,  or  matched  with  other  similar  numeric  simulations  to  help  establish  the 
accuracy  of  the  code.  Solutions  to  Poisson’s  equation  for  the  potential,  electron  and  ion 
continuity  equations,  and  the  electron  energy  equation  are  solved  semi-implicitly  in  a 
sequential  manner.  Each  of  the  governing  equations  is  solved  by  casting  them  onto  a 
tridiagonal  grid,  and  using  the  computationally  efficient  Thomas  algorithm  to  solve  ID 
regions  in  a  single  iteration.  The  Scharfetter-Gummel  flux  discretization  method  is  used 
to  add  stability  to  the  code  when  transitioning  from  a  field  to  diffusion  dominated  region 
or  vice  versa.  Estimates  for  the  ionization  and  recombination  rates  and  for  the  transport 
coefficients  of  the  background  gas  are  calculated  as  a  function  of  the  local  average 
electron  energy,  and  updated  for  every  calculation  point  in  the  domain  on  the  completion 
of  the  solution  to  the  electron  energy  equation.  Results  are  then  recorded  and  used  as  a 
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starting  point  for  the  next  time  step.  The  basis  for  the  2D  geometry  relies  heavily  on  its 
ID  counterpart,  but  with  the  inclusion  of  a  Gauss-Seidel  line  iterative  method  to  solve 
Poisson’s  equation,  and  a  superposition  method  for  keeping  track  of  ion,  electron,  and 
electron  energy  fluxes.  Appropriate  boundary  conditions  are  implemented  to  close  the 
computational  region  at  the  boundaries,  take  into  account  any  space  charge  present  in  the 
computational  area,  and  keep  track  of  surface  charge  buildup  on  any  dielectric  surfaces  if 
present.  In  the  2D  geometry,  momentum  transfer  to  the  background  gas  could  be 
estimated  by  tracking  the  charge  particle  movement  and  their  interactions  with  the  neutral 
background  gas;  and  the  effects  of  various  electrode  placements,  driving  potentials,  and 
driving  frequencies  (if  an  AC  potential  source  is  present)  could  be  used  to  further  the 
understanding  of  DBD  operation.  Results  from  the  ID  and  2D  code  are  presented,  as 
well  as  limitations  of  the  current  work  and  proposed  future  refinements. 
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COMPUTATIONAL  MODELING  OF  THE  DIELECTRIC  BARRIER 
DISCHARGE  (DBD)  DEVICE  FOR  AERONAUTICAL  APPLICATIONS 


I.  Introduction 


Background  and  Importance  of  the  Dielectric  Barrier  Discharge 

The  Dielectric  Barrier  Discharge  (DBD)  device  has  been  put  to  use  since  1857 
when  Werner  von  Siemens  used  to  produce  ozone  at  atmospheric  pressures  from  the 
oxygen  in  air.  Today,  DBD  devices  are  still  in  use  as  ozone  generators,  but  have  found 
many  more  uses  in  industry,  science,  and  military  applications.  Such  uses  include  plasma 
television  displays,  pumps  for  C02  lasers,  toxic  gas  decomposition  for  pollution  control, 
and  the  modification  of  air  flow  over  a  surface  (14:1819).  The  Fig  1  shows  the  nominal 
flow  control  set-up  for  the  DBD  device  as  described  by  Font  et  al  (8: 1). 


alternating  voltage 


Plasma  actuator  configuration. 


Figure  1.  DBD  device  configuration  for  aeronautical  application  (8:1). 
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This  technology  is  of  particular  interest  to  the  United  States  military,  specifically 
the  Air  Force,  to  improve  the  effectiveness  of  their  aircraft  and  weapon  systems.  DBD 
devices  have  been  shown  to  induce  air  flow  over  a  surface,  reattach  flow  to  an  airfoil  at 
high  angles  of  attack,  such  as  shown  in  Fig  2;  and  show  promise  of  eventually  replacing 
flaps  and  ailerons  on  aircraft  (DBD  devices  are  commonly  called  plasma  actuators  when 
referring  to  flow  control).  By  eliminating  the  hydraulics  on  aircraft  and  replacing  them 
with  electrical  wires,  plasma  actuators  would  open  the  door  to  building  stronger,  lighter, 
more  robust  airfoils,  with  faster  maneuvering  response  time  and  fewer  moving  parts  to 
malfunction.  This  is  especially  advantageous  when  considering  Unmanned  Combat 
Aerial  Vehicles  (UCAVs),  and  pushing  the  performance  envelope  above  and  beyond 
what  human  physiology  can  withstand. 


Plasma  off 


Plasma  on 


Figure  2.  Boundary  layer  flow  re-attachment  caused  by  a  DBD  device  (17:2124). 


The  maturation  of  DBD  actuator  technology  by  the  United  States  represents  the 
next  generation  of  design  enhancements  that  will  continue  to  keep  Air  Force  jets  flying 
higher,  faster,  and  farther  with  the  potential  of  being  stronger,  lighter,  and  more  highly 
maneuverable.  The  development  of  robust  and  accurate  computational  code  for 
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simulating  the  outcome  of  these  ongoing  experiments  would  represent  a  great  step  toward 
understanding  their  behavior  more  completely.  Replacing  the  experimental 
configurations  that  need  to  be  built  and  tested  in  a  laboratory  with  a  numeric  simulation, 
would  result  in  great  saving  in  time  and  resources,  and  could  possibly  lead  to  faster 
development.  In  short,  an  effective  and  efficient  numerical  code  would  help  mature  and 
deploy  the  DBD  actuator  technology  more  quickly  than  traditional  methods  and  also 
result  in  a  significant  costs  savings  to  the  Air  Force. 

The  focus  of  this  thesis  will  be  to  provide  a  greater  understanding  of  the  discharge 
dynamics,  through  observing  momentum  and  energy  transfer  in  the  plasma  using  fluid 
dynamics  and  computational  methods.  Eventually,  the  availability  of  accurate  modeling 
software  that  is  grounded  in  theory  will  allow  engineers,  once  left  with  only  trial  and 
error,  to  optimize  their  designs  in  a  cost  effective  way.  A  robust  and  accurate 
computational  tool  would  help  predict  the  energy  and  momentum  imparted  to  the 
surrounding  airflow  via  the  DBD  actuator  for  different  device  configurations  and 
operating  parameters.  This  will  provide  the  Air  Force  with  a  low  cost,  quick  turn  around, 
effective  way  of  designing  their  next  generation  of  aerial  vehicles. 

DBD  Device  Configuration 

A  DBD  device  can  be  made  from  any  configuration  of  electrodes,  separated  by  a 
dielectric  barrier,  operated  under  an  alternating  current  (AC)  configuration,  as  opposed  to 
direct  current  (DC).  The  devices  are  commonly  used  in  the  glow  discharge  region  of  the 
plasma  spectrum,  where  the  number  density  of  positively  and  negatively  charge  species 
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are  roughly  equivalent,  and  a  neutral  gas  many  orders  of  magnitude  greater  in  density  is 
present.  This  background  pressure  can  range  from  a  few  Torr  to  1  Atmosphere  (ATM) 
(20:8).  Initially,  the  applied  potential  difference  on  the  electrodes  must  be  great  enough 
to  initiate  gas  breakdown.  While  operation  is  possible  for  a  short  time  using  a  DC 
potential,  the  discharge  will  eventually  extinguish  due  to  a  buildup  of  charge  on  the 
dielectric  barrier  as  notionally  seen  in  Fig  3.  This  buildup  effectively  reduces  the  applied 
potential  across  the  gap  between  the  electrodes  and  the  discharge  is  extinguished.  This 
buildup  may  play  a  crucial  role  in  the  continued  operation  and  characterization  of  these 
devices  under  AC  operating  conditions.  It  might  act  to  seed  the  ionization  of  the  neutral 
gas  on  the  reverse  cycle  of  the  potential  as  shown  below  in  Fig  3.  This  seeding  could  in 
fact  be  the  mechanism  by  which  an  asymmetric  force  on  the  neutral  flow  is  coupled  to  the 
background  gas  during  operation  (9:9).  For  this  project,  the  charge  will  be  considered 
trapped  on  the  dielectric  surface  and  only  secondary  emission  from  exposed  electrode 
surfaces  due  to  ion  impact  will  be  characterized. 


Figure  3.  Charge  build-up  on  a  DBD  dielectric  under  AC  operating  conditions  (7:2). 
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The  dielectric  barrier  can  be  either  an  insulating  coating  on  the  anode  or  cathode 
of  the  device,  or  a  layer  of  insulating  material  placed  somewhere  in  between.  A  variety 
of  DBD  configurations  exist,  some  of  which  are  shown  in  Fig  4.  These  various 
arrangements  stem  from  the  multiple  uses  of  DBD  type  devices  for  industrial 
applications,  plasma  T.V.  display  light  generation,  and  aeronautical  flow  control 
purposes.  It  should  be  noted  that  while  these  DBD  devices  have  found  widespread  usage 
in  the  modern  world  because  of  their  proven  operational  capability,  the  physics  behind 
their  operation  and  plasma  dynamics  involved  are  not  well  understood.  Therefore,  this 
thesis  will  focus  on  characterizing  the  mechanisms  involved  in  DBD  operations  through 
the  use  of  macroscopic  fluid  equations.  The  numeric  simulation  will  be  done  for  both  ID 
and  2D  geometries. 


High 
Voltage 
AC  I 
Generator 


High  Voltage 
Electrode 


Dielectric 

Barrier 


Discharge 


Ground 

Electrode 


Figure  4.  Various  geometrical  arrangements  for  the  DBD  device  (14: 1820). 
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II.  Computational  Modeling 


Introduction 

Tracking  the  position  of  all  the  individual  particles  in  the  DBD  system,  as  well  as 
their  ancillary  velocity  vectors  for  any  system  with  more  than  a  few  constituents  is 
problematic,  if  not  impossible.  Of  course,  the  only  exact  solution  is  to  track  each 
constituent  on  a  microscopic  scale;  however,  memory  requirements,  as  well  as 
computational  time  limitations  force  an  approximate  solution  where  only  the  macroscopic 
parameters  of  the  system  are  accounted  for  and  varied.  In  making  a  few  physically 
reasonable  approximations,  the  fluid  equations  governing  the  DBD  system  can  be  derived 
from  the  Boltzmann  equation,  allowing  for  numeric  modeling  of  different  conditions  and 
configurations  in  a  reasonable  amount  of  time  without  losing  too  much  information  about 
the  physical  reality  of  the  model.  For  this  thesis,  macroscopic  fluid  type  equations  will 
be  used  to  track  the  location  and  movement  of  the  electrons  and  positively  charged  ions 
in  the  system,  as  well  as  the  average  local  electron  energy  density  at  each  spatially 
discretized  location. 


Moments  of  the  Boltzmann  Equation 

The  Boltzmann  equation  is  defined  by: 


B 


df_ 

dt 


+  v  •  V/  +a-V-f  =  I 

v  1  St 


/  collision 


(1) 


This  equation  represents  the  temporal  evolution  of  the  distribution  function,  f  (or  the 
probability  of  finding  a  plasma  species  with  a  particular  position  and  velocity  ( v  )).  The 
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first  term  represents  the  change  in  the  distribution  function  over  time;  this  is  ultimately 
what  will  be  solved  for  to  give  a  solution  over  all  time  to  a  system  given  some  set  of 
initial  conditions.  The  second  term  describes  how  the  distribution  changes  in  space,  or 
configuration  coordinates,  as  the  particles  move  through  the  system  with  a  certain 
velocity  and  direction.  The  third  shows  how  the  species  distribution  is  changing  due  to 
external  forces  acting  to  accelerate  the  charged  particles.  And  finally,  the  tenn  on  the 
right  hand  side  represents  collisional  interactions  among  the  plasma  species  as  time 
passes,  such  as  ionization  and  recombination,  and  how  these  affect  the  total  distribution 
of  the  system.  The  zeroeth  and  second  moments  of  the  Boltzmann  equation,  are  defined 
by: 


M o m  ent0 


oltzmann 


dv 


(2) 


and 


M  o  merit  2 


oltzmann 


—  m  svv  dv 
2 


(3) 


respectively.  The  integrals  are  done  over  all  space  with  respect  to  the  velocity 
coordinates.  These  two  moments,  along  with  the  solution  to  the  first  moment,  also  called 
the  momentum  equation,  which  will  be  approximated  using  the  Scharfetter-Gummel  flux 
discretization  technique  described  later,  can  be  solved  to  yield  the  continuity  and  energy 
equations  for  a  particular  species,  with  m  designating  the  mass  of  the  species.  The 
species  (designated  by  a  superscript  s)  included  in  this  thesis  will  be  limited  to  singly 
positively  charged  ions  and  electrons  only.  A  brief  development  of  the  continuity 
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equation  will  be  looked  at  next,  but  results  for  the  energy  equation  will  only  be  stated. 
Full  derivations  can  be  found  in  almost  any  introductory  plasma  physics  course  text 
(6:19-26). 

By  simplifying  each  tenn  of  the  zeroeth  moment  of  the  Boltzmann  equation 
defined  above,  the  continuity  equation  can  be  shown  to  develop  as  follows  (19): 


J  f +v.V/  +  «.V,/ 


Jv-V/rfv  +ja-'V~fdv  =| 


collision 


Term  1  Term  2  Term  3  Term  4 


Term  1  reduces  to: 


dn 

dt 


because  n  =  /  dv ,  is  the  number  density  of  the  species. 


Term  2  reduces  to: 


(5) 


dv  V  •  v  =  V  •  J  /  dv  v  =  V  •  nv  =  V  f  (6) 

recognizing  that  the  number  density  times  the  velocity  is  defined  as  the  species  flux  (T  ) 
into  or  out  of  any  particular  volume  element. 

Term  3  reduces  to: 
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(V) 


ja  ■  V-/  dv  =  |V-  •  a  f  dv  =  |  af-ndS-0 

Surface 

by  using  the  divergence  theorem,  and  recognizing  that  any  distribution  goes  to  zero  as  the 
surface  used  in  the  integral  approaches  infinity. 

Finally,  term  4  reduces  to: 


resulting  in  a  source  and  loss  tenn  detennined  by  the  collisional  interaction  of  the 
particles  in  the  system. 

The  resulting  continuity  equation  can  be  rewritten  from  Eq  (4)  as: 

— +  V-f  =S-L.  (9) 

dt 

If  ionization  and  recombination  of  the  positively  charged  ions  and  electrons  are 
the  only  source  and  loss  terms  considered  for  this  project,  then  the  S  —  L  tenn  above  can 
be  rewritten  as: 

S  -  L  =  vnc  -  j3nen‘  (10) 

with  v  representing  the  electron  impact  ionization  frequency  on  the  neutral  gas,  and  the 
recombination  coefficient  represented  by  p.  The  local  temporal  change  of  the  number 
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density  and  the  divergence  of  the  flux  in  ID,  along  the  X-axis  in  this  case,  can  be  recast 
in  a  finite  difference  fonn  as: 


'■N  S.t  S,t— 1 

on  _n  —n 
dt  At 


(11) 


and 


V-f 


(12) 


respectively.  The  superscript  t  refers  to  a  specific  point  in  time.  The  divergence  term 
reduction  above  is  commonly  referred  to  as  a  three-point  formula,  even  though  the  third 
point  does  not  appear  in  the  equation  (4:171-172).  This  ID  discretization  is  appropriate 
for  this  project  since  particle  fluxes  into  or  out  of  a  local  volume  will  only  be  considered 
along  one  specific  axis  at  a  time.  Note  that  the  flux  characters  retain  their  vector  quality, 
even  when  reduced  to  ID  in  the  example  above;  this  indicates  the  ability  of  the  flux  to  be 
directed  in  either  a  positive  or  negative  sense. 

Finally,  all  tenns  from  the  previous  development  can  be  combined  to  yield  the  ID 
continuity  equation  for  this  project: 


n 


_  S  ,t — 1 


At 


■  + 


right 


_ 1  left 


Ax 


vn 


e,t- 1 


—  fin 


e,t- 1  i.t—1 

n 


(13) 


The  electron  energy  equation  takes  on  a  similar  form: 
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(14) 


e,t  e,t 

n  u 


e,t— 1  e,t- 1 

—  n  u 


At 


5  r 

-+ — 

3 


right 


pfl  u 
~  1  left 


Ax 


=  -qTe 


■E-NkLn 


e,t- 1 


neue  designates  the  electron  energy  density.  The  first  tenn  on  the  right  hand  side  is  a 
source  due  to  Joule  heating,  with  q  as  the  unsigned  charge  on  an  electron;  and  the  second 
tenn  is  a  collisional  power  loss  associated  with  electron-neutral  collisions.  If  the  units  of 
energy  are  specified  to  be  in  eV  instead  of  the  S.I.  unit  of  Joules,  then  the  q  in  Eq  (14) 
above  is  divided  out  in  the  conversion,  and  everything  else  remains  the  same.  Neutral 
number  density  is  determined  by  specifying  the  pressure  ( P )  at  which  the  simulation  will 
be  running,  and  converting  via: 


N  = 


(15) 


with  T  as  the  temperature  of  the  system  and  kb  as  Boltzmann’s  constant. 

Discretization  of  Equations  and  the  Scharfetter-Gummel  technique 

With  the  fluid  equations  specified,  the  next  step  is  to  approximately  describe  them 
for  the  discrete  points  in  space  that  make  up  the  computational  domain.  In  a  2D 
geometry,  but  considering  only  the  change  in  density  due  to  movement  of  the  species 
along  the  X-axis  in  this  example,  the  electron  and  ion  continuity  equations  can  be  written 
as: 


sJ  s,/-l 

—  n 

x,y  x,y 


At 


1  x+Yi^y  L  x~Y>y 


Ax 


vn 


-fin 


n 

x,y  x,y 


(16) 
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While  the  electron  energy  equation  can  be  written  as: 


e.t  _  e.t 

n  u 

x,y  x,y 


x,y  x,y 


At 


+  ■ 


sr: 


-r 


/2’y  _ 


Ax 


retE 


e,t— 1 


X,y 


NkLnx,y 


(17) 


with  the  subscripts  now  designating  a  particular  position  in  2D  space.  This  project  will 
set  the  distance  between  evaluation  points  to  be  equal,  that  is  Ax  =  Ay  .  The  final  form  of 
these  equations  suggests  that  for  given  initial  conditions  the  evolution  of  the  charged 
species  densities  and  energy  density  can  be  solved  for  any  later  time  in  a  sequential 
manner.  Solving  then  for  the  current  time  step,  these  two  equations  can  be  rewritten  as: 


n. 


x,y 


-(S'1 -i: 


*,y 


T~ 'S,t  T" 'S,t 

x+y2,y  x~Yi,y 

Ax 


)At  +  n 


5,^-1 

x,y 


(18) 


and 


e.t  „  e.t 

nu 

x,y  x,y 


=  (S' 


■L” 


5  r 


x+y,y 


-r 


x,y 


Ax 


)At  +  n 


e,t- 1  e,t- 1 

u 

x,y  x,y  • 


(19) 


Having  described  the  source  and  loss  terms  earlier,  the  flux  terms  in  the  above 
equations  will  now  be  described  using  the  Scharfetter-Gummel  flux  discretization 
scheme.  This  flux,  which  appears  in  both  the  continuity  and  energy  equations,  is 
obtained  from  the  solution  of  the  first  moment  of  Boltzmann’s  equation  by  applying  the 
drift  and  diffusion  approximation  (12:3).  By  using  this  flux  description  for  the  continuity 
and  energy  equations,  computational  stability  is  added  to  the  model,  allowing  for 
relatively  coarser  grid  point  spacing.  This  is  a  result  of  the  ability  of  the  Scharfetter- 
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Gummel  technique  to  smooth  the  transition  between  field  driven  and  diffusion  dominated 
regions  (12:4).  This  smoothness  is  achieved  by  the  built-in  ability  of  the  flux 
representation  to  account  for  the  direction  of  the  electric  field  at  any  given  point  {up- 
winding),  relate  information  about  the  species  densities  on  either  side  of  that  point,  and 
appropriately  weight  the  contributions  of  each  term  when  describing  the  flux  into  or  out 
of  the  point.  The  up-winding  effect  of  the  Scharfetter-Gummel  technique  is  not  easily 
seen  at  first,  and  a  more  in-depth  example  describing  the  behavior  under  limiting 
conditions  can  be  found  in  the  work  by  Hilbun  (13:9-10).  The  flux  terms  for  the 
continuity  equations  are: 


•CDS'  \ 

^  ys,t  ^ 

Ar  j 

texp(Z^-Dj 

(20) 


and 


f<-i,,^exp(Z^ 

)-  ns''Ds''-x\ 

,y  /  x,y  x,y 

(  ZW  1 

x-'A-y 

v  Ax 

J 

exp (Zs\,  -1) 

V  ^  x~A-y  J 

(21) 


Where  the  Z  terms  are  defined  as: 


=  ~sign[s]  ^v+,_r  (# 
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S,t- 1 

x+y2,y 


x+l,y 


(22) 
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~sign[s] 
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.s,t- 1 

_ 'A± 

DSJ 
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(&.V-&- 1.,) 
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(23) 
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with  sign[s]  taking  the  value  of  +1  for  the  positive  ions  and  -1  for  electrons.  The  (f)  ’s 


represent  the  scalar  potential  at  a  given  grid  point,  and  D  and  //  are  the  diffusion  and 
mobility  coefficients  respectively.  Special  limiting  if  statements  were  included  in  the 
computational  code,  to  account  for  regions  where  the  change  in  potential  was  very  small 
or  zero,  or  where  the  change  was  very  large.  This  kept  the  F  terms  defined  below  from 
returning  as  undefined,  or  as  a  number  greater  than  the  computer  could  handle  due  to  the 
exponential  factors  in  the  terms  of  Eq  (24).  It  should  be  noted  that  the  values  of  Z,  D,  and 
ju  for  the  electron  energy  flux  are  the  same  as  those  appearing  in  the  electron  continuity 

equation  (i.e.  z£'-'=Z£1,  =  Df;\  /C’^=/0-  To  further  reduce  the 


equations  into  a  more  manageable  form,  the  FI  and  F2  terms: 


F\[Z]  = 


F2[Z\  = 


Z  exp  (Z) 
exp(Z) - 1 
Z 

exp(Z) - 1 


(24) 


are  introduced.  Shown  in  the  figure  below  is  a  plot  of  these  two  F  functions. 


Figure  5.  Plot  of  FI  and  F2  as  functions  of  Z. 


14 


The  limit  of  the  F  functions  for  small  Z’s  in  the  range  of  +1CT7  about  zero  was  set 
to  1  based  on  L’Hopital’s  rule,  while  the  limit  for  the  Z’s  that  were  greater  than  150  or 
less  than  -150  was  set  equal  to  Z  and  0  for  FI  and  F2  respectively.  After  a  little  algebra, 
new  forms  for  the  flux  terms  can  be  shown  reduced  as: 


yD*,y 


lFl[Zs;‘]-K‘D 


s,t- 1 


F2[Z:tyJ 


Ax 


(25) 


and 


OP4] 


5,?-l 


*x-l,y^x-l,y 


ns’  D 

x,y  x,y 


F2  [z;^y\ 


Ax 


(26) 


Finally,  substitutions  for  the  previously  determined  representations  for  source  and 
loss  can  be  made.  The  final  form  of  the  electron  continuity  equations  in  ID  along  the  X- 
axis  reduces  to: 


<y  ((A*)2  +  D«;'F\[Z:i/  y  ]M  +  D»~'F2[Z:yr  ]  A/)  + 

<:y~DXx2[zyt!,]&t)  = 


1 


1 


(Ax)  +  -  vx  ynl;  ~  At(Ax)  -  -  pnl;~  Af(Ax)" 


(27) 


The  ion  continuity  equations  in  ID  along  the  X-axis  can  be  re-written  as: 
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(28) 
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where  the  source  and  loss  terms  for  this  project  will  always  depend  on  the  previous 
values  of  the  electron  and  ion  densities  to  ensure  generation  or  recombination  of  the 
charged  species  occurs  in  pairs,  making  the  numeric  code  is  stable,  and  not  allowing  any 
inequality  between  production  and  loss.  In  a  similar  manner,  the  electron  energy 
equation  in  ID  along  the  X-axis  takes  on  the  final  fonn: 
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(<A*)! ) - (f  1„  •  £,„)  A/( Axf  -  i n«-'kLNAt( Ax? 


(29) 


where  the  electron  flux  term,  P  ,  at  the  grid  point  is  the  average  of  the  electron  flux  at 

the  current  time  step  of  the  closest  mid-points  which  is  solved  for  during  the  solution  of 
the  electron  continuity  equation  using  the  Scharfetter-Gummel  discretization  technique. 
These  flux  values  are  determined  when  the  electron  continuity  equation  is  solved.  It 
should  be  recognized  that  the  values  of  the  electron  flux  and  electric  field  in  the  energy 
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equation  retain  their  vector  nature,  indicating  that  they  lie  along  a  particular  axis  in  2D 
space.  The  factor  of  14  multiplying  the  source  or  loss  terms  in  all  the  governing  equations 
is  a  manifestation  of  calculating  the  source  or  loss  for  each  volume  element  in  2D  space 
as  a  superposition  of  both  ID  elements  along  the  two  possible  axes.  This  will  be 
elaborated  on  more  in  Chapter  3. 

With  the  three  governing  equations  now  reduced  into  a  compact  tridiagonal  form, 
the  evolution  of  the  ion  and  electron  densities  and  the  electron  energy  from  some  initial 
condition  can  now  be  found,  given  that  appropriate  values  for  the  boundary  conditions  of 
the  system  and  the  existing  unknowns  can  be  detennined. 

Boundary  Conditions 

Since  the  Scharfetter-Gummel  flux  representation  depends  on  neighboring  points 
to  detennine  intermediate  values,  and  this  model  will  assume  that  the  space  charge 
density  at  any  surface  boundaries  is  zero,  realistic  conditions  for  the  fluxes  at  these 
boundaries  need  to  be  defined  (2:1379). 

The  flux  of  the  electrons  at  a  surface  boundary  will  be  defined  as: 

Y'eJ  _ |_  1  e,t  ethennal.t  —  . 

boundary  ^  adjacent  adjacent  T  /  adjacent  (30) 

where  the  plus  or  minus  terms  indicates  that  the  thermal  flux  tenn  is  always  directed 
towards  the  surface  boundary,  and  that  the  secondary  emission  is  directed  away  from  the 
boundary,  y  represents  the  secondary  emission  coefficient  for  positive  ion  impact  onto  the 
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surface  of  any  exposed  electrode,  y  will  be  set  to  0.02  for  this  project.  is  the 

value  of  the  thermal  velocity  of  the  electrons  at  the  point  next  to  the  boundary,  expressed 
in  terms  of  the  average  electron  energy  at  that  point: 


e  thermal 
^ adjacent 


,t 


4.19x10s 


e,t- 1 

adjacent  * 


(31) 


The  average  electron  energy,  ue,  is  found  by  dividing  the  energy  density  in  the  solution  of 
the  energy  equation,  by  the  number  density  in  the  solution  of  the  electron  continuity 
equation. 

The  conditions  for  the  ion  flux  at  the  surface  boundaries  are  somewhat  different, 
but  still  give  a  realistic  value  over  the  time  scale  used  because  of  their  much  slower 
response  to  any  applied  fields  and  relatively  cooler  temperature  when  compared  to  the 
electrons.  That  being  said,  the  movement  of  the  ions  is  modeled  as  field  driven  only  at 
the  surface  boundaries  using  the  Scharfetter-Gummel  technique.  To  satisfy  the  condition 
of  having  only  field  driven  flux,  the  diffusion  dependency  of  the  technique  must  be 
zeroed.  This  is  done  by  temporarily  setting  the  ion  density  at  the  surface  boundary  equal 
to  its  adjacent  grid  point.  One  final  if  statement  must  be  incorporated  however.  The 
boundary  flux  must  only  exist  if  the  electric  field  is  also  directed  toward  the  boundary. 
This  condition  keeps  ions  from  moving  out  of  the  region  that  was  temporarily  assigned  a 
value  to  zero  the  diffusion  term,  but  is  otherwise  modeled  as  having  no  charge  density. 
The  final  form  of  the  ion  flux  at  the  surface  boundary,  which  can  be  integrated  into  the 
previously  detennined  ion  continuity  equations,  is: 
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_  (F\[Ziidjacent±y2  ] 

adjacent 

where  the  ±  term  in  the  Z  factor  indicates  if  an  adjacent  grid  point  is  either  to  the  left  or 
the  right  of  the  surface  boundary  in  this  X-axis  example. 

Finally,  the  electron  energy  flux  into  the  surface  boundary  takes  on  a  form  very 
similar  to  that  of  the  electron  flux.  It  is  modeled  as: 


^  adjacent'/,  ])n 


i,t 
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Any  electron  energy  flux  out  of  the  boundary,  due  to  secondary  emission,  is  accounted 
for  in  the  Joule  heating  term  of  the  energy  equation. 

Conditions  for  non-surface  boundaries,  for  the  2D  geometry,  are  implemented  by 
using  ghost  cells  around  the  computational  domain  that  have  the  same  value  for  densities 
as  their  neighboring  cells,  and  where  the  perpendicular  component  of  the  electric  field 
has  been  defined  to  be  zero,  so  as  not  to  allow  any  loss  in  the  system.  While  the  zeroing 
of  the  perpendicular  component  of  the  electric  field  at  the  non-surface  boundary  is  not  a 
physical  reality,  it  is  a  good  approximation  if  the  size  of  the  domain  and  distance  of  the 
boundary  is  far  enough  away  from  the  active  region  of  the  simulation  (1 1:40). 

BOLSIG  Fits  for  Rate  and  Transport  Coefficients 

With  the  boundary  conditions  adequately  defined,  the  values  of  the  transport 
parameters  and  kinetic  rate  coefficients  in  the  continuity  and  energy  equations  will  now 
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be  defined.  The  values  of  /?,  y,  and  ju'/y  or  D'/  are  assumed  to  be  constants  that  were 

measured  experimentally  given  the  material  make-up  of  the  system,  the  type  of  gas  in  the 
system,  and  the  pressure  of  the  gas  in  the  system.  Using  the  Einstein  relation  to  define  a 
relationship  between  ju'/  and  D['\.  ,  for  the  positively  charged  ions  in  the  system,  it  can 
be  shown  that: 


Dl 


x,y 


M: 


l,t 

'x,y 


khT 


q 


(34) 


Approximating  room  temperature  for  this  simulation  leaves  a  value  of  about  1/40 
for  the  right  hand  side  of  the  equation,  showing  that  if  one  parameter  value  is  known 
from  an  experimental  result,  the  other  can  be  easily  detennined. 

The  values  of  ,  De/ ,  vx  Y ,  and  kL  as  a  function  of  local  average  electron  energy 

can  be  determined  by  using  the  BOLSIG  Boltzmann  equation  solver.  This  excellent 
freeware  tool  has  proven  reliable  and  accurate,  as  well  as  easy  to  use.  Curve  fits  of  the 
unknown  coefficients  as  a  function  of  average  electron  energy  can  be  obtained  for 
specified  gases  and  mixtures  of  gases,  as  well  as  their  dependence  on  ambient  pressure. 
This  project  uses  the  BOLSIG  results  for  a  pure  atmosphere  of  either  Ar  or  N2.  Two 
representative  plots  of  the  numerical  results  from  BOLSIG,  as  well  as  their  associated  fits 
are  shown  below  in  Fig  6  and  Fig  7  (3). 
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N2  BOLSIG  Data  -  Mobility  x  N  vs.  Energy 


Average  Electron  Energy  (in  eV) 


Figure  6.  Computational  fit  for  the  electron  mobility  as  a  function  of  average  energy. 


N2  BOLSIG  Data  -  Loss  Coeff  vs.  Energy 


Average  Electron  Energy  (in  eV) 


Figure  7.  Computational  lit  for  the  kL  coefficient  as  a  function  of  average  energy. 


In  the  first  graph,  the  fit  obtained  from  BOLSIG  was  determined  to  be  accurate 
enough  to  use  for  this  project,  the  second  was  fitted  to  a  more  appropriate  polynomial  to 
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more  accurately  simulate  the  discreet  BOLSIG  results.  These  fits  were  calculated  for  the 
range  0. 0-5.0  eV  and  were  determined  to  be  in  good  agreement  with  the  calculated  values 
in  the  range  0.3-4.75  eV,  therefore  upper  and  lower  boundaries  were  placed  on  the 
average  electron  energy,  as  to  not  incorrectly  describe  the  rate  and  transport  coefficients 
by  using  the  fits  to  calculate  them  in  a  region  were  they  are  not  appropriately  represented. 
An  example  of  this  limiting  factor  can  be  seen  in  the  behavior  of  the  fit  for  the  energy 
loss  coefficient  above  where  the  value  of  the  polynomial  fit  drops  below  zero  at  around 
0.3  eV. 


The  last  set  of  unknowns,  potential  and  electric  field  values,  for  all  points  on  the 
computational  grid  must  now  be  detennined  in  order  to  arrive  at  a  solution  of  the 
continuity  and  electron  energy  equations.  Since  the  computational  set-up  will  include 
regions  of  different  permittivity,  Gauss’s  law  in  the  presence  of  dielectrics  must  be  used 
to  accurately  account  for  any  bound  charge  which  might  contribute  to  the  polarization  of 
the  material  under  an  applied  electric  field  (10:175-180).  This  law  in  equations  form  is: 

V-D  =  pf.  (35) 

D  is  the  electric  displacement  vector  defined  by  D  =  s0E  +  P ,  where  £0  is  the 

pennittivity  of  free  space,  E  is  the  electric  field  in  the  material,  and  P  is  the  dipole 

moment  per  unit  volume.  The  vector  P  is  used  to  correctly  account  for  the  bound  charge 
density  associated  with  polarization  by  an  externally  applied  electric  field  in  a  dielectric 

material.  In  linear  media  P ,  can  be  further  defined  as  Ps  s0  %eE ,  allowing  for  a  new 
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equivalent  definition  of  the  displacement  vector:  D  =  e0(l  +  %e)E  .  Finally,  making  the 
substitution  s  =  s0( \  +  %e)  Gauss’s  law  can  be  rewritten  as: 


V-sE  =  pf  (36) 

where  s  is  the  pennittivity  of  the  material.  Defining  this  pennittivity  as  a  multiple  of  a 
relative  permittivity  and  the  pennittivity  of  free  space  gives: 

£  =  £reo  •  (37) 


This  relative  permittivity,  sr ,  also  referred  to  as  the  relative  dielectric  constant  of  the 

material,  is  well  known  for  many  materials  through  experimental  measurement.  One  of 
the  most  common  dielectric  materials  used  in  the  aerospace  industry  today  is  that  of 
Kapton®  tape  manufactured  by  Dupont,  with  a  relative  dielectric  constant  of  4.0;  and  will 
therefore  be  the  relative  dielectric  constant  chosen  to  define  the  permittivity  of  the 
dielectric  region(s)  between  the  anode  and  cathode  in  this  project. 

Using  the  three-point  formula  to  solve  for  the  divergence  of  the  electric  field  in 
Eq.  (36)  above,  and  recalling  that  the  dielectric  material  being  used  exhibits  linear 
properties,  Gauss’s  law  can  be  rewritten  in  2D  as: 
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(38) 
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The  free  charge  on  the  right  hand  side  of  the  equation  is:  p*  =  q(nlf  1  -  nx‘y  1 ) .  The 

electric  field  component  in  Eq  (38)  above  can  be  defined  in  terms  of  the  negative  gradient 
of  the  potential  as: 

E  =  -  V (j) .  (39) 


Therefore,  substitutions  can  be  made  for  the  individual  components  in  Gauss’s  law  as: 
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Furthermore,  using  the  previously  defined  assumption  for  this  project  that  Ax  =  Ay ,  the 
equation  reduces  to: 

“V/24h  +(%>'  +Sx-^y  +  V/2  +V/2Mv  -£x,^Ay+\  =Px,JM  +£x$jL, y  +8x-%y<i>xA,y  -(41) 


Combining  like  terms,  and  assuming  that  the  values  of  s  are  equal  over  the  region  of 
interest,  leaves  the  final  fonn  of  Gauss’s  law  as: 


px  (Ax)2 

x,y- 1  +  40x,y  ~  0x,y+l  =  ~ ~  +  <4+1, y  +  <4-1 ,v  ' 


(42) 
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However,  if  the  values  of  s  vary  over  the  spatially  discrete  grid,  meaning  a 


dielectric  interface  is  encountered,  Gauss’s  law  take  on  the  form: 
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This  is  because  the  value  of  pf  at  any  of  these  points  for  this  project  is  defined  as  zero, 
the  value  of  e  at  the  dielectric  interface  has  been  defined  as: 


^  =  <44) 

Any  accumulated  free  charge  on  the  dielectric  surface  has  been  accounted  for  by 
substituting  ax  =  px  Ax  into  Gauss’s  law.  Modeling  the  dielectric  in  the  system  as 

sticky,  that  is  any  flux  of  plasma  species  incident  on  the  surface  adheres  and  remains  at 
the  location  of  incidence  and  assuming  that  any  neutralization  of  the  charge  species  by 
recombination  takes  place  instantly  at  the  surface,  a  can  be  defined  as  the  accumulation 
over  time  of  any  fluxes  incident  on  the  surface.  The  surface  charge  density  will  be 
calculated  using: 


_t  t- 1  ,  / 

=Crx  +</[ 


1 


(45) 


reflecting  the  fact  that  the  dielectric  surface  for  the  2D  geometry  in  this  project  will  be 
constrained  to  lie  exclusively  along  the  X-axis.  A  more  accurate  approximation  of  the 
surface  charge  accumulation  and  relaxation,  could  be  obtained  by  adding  the  term: 
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(  ^ conductivity^ interface  '  (46) 

to  the  previous  equation  if  the  charge  on  the  surface  was  negative,  which  represents  the 
rate  at  which  electrons  will  flow  through  the  dielectric,  or  the  dielectric  relaxation  time. 
However,  if  this  term  was  included,  a  more  complicated  description  of  the  surface  charge, 
including  emission  from  the  dielectric  under  appropriate  field  conditions,  would  have  had 
to  be  considered.  Therefore,  the  choice  was  made  to  simply  consider  any  charge  incident 
on  the  surface  as  sticky  and  not  allowed  to  leave,  only  recombine. 

The  previously  defined  equations  for  the  potential  at  individual  grid  points  can  be 
united  to  solve  an  entire  line  of  potential  values  in  space  simultaneously  by  casting  them 
into  a  matrix  of  tri diagonal  form.  An  example  of  this  matrix  for  a  ID  computational 
region  divided  into  N  =  9  potential  points,  with  potential  points  0  and  8  being  specified 
and  one  dielectric  interface  at  point  4  is  shown  below. 

(2-io  o  oo  0Y4A  (  ^  ^ 

-1  2  -1  0  0  0  0  (j)xl  0 

0  -1  2  -1  0  0  0  (f>xi  0 

0  0  ~Sx,y-/  £^/+i W2  “Yr+h  0  0  4,4  =  AA  (47) 

0  0  0  -1  2  -1  0  4,5  Pjte? 

0  0  0  0  -1  2  -1  4,6  A,6(^)2 

v0  0  0  0  0  -1  2  JU,7 J  +4.8, 

The  equation  above  corresponds  to  the  ID  set-up  shown  below  in  Fig  8.  Squares  denote 
locations  at  which  the  potential  is  calculated  as  well  as  the  value  of  any  free  charge  in  the 
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system.  Surface  charge  is  allowed  to  accumulate  on  the  dielectric  surface  for  this 
example.  X’s  indicate  where  the  fluxes  between  points  will  be  calculated. 


ID  Dielectric  Barrier  Discharge  Set-up  for  N=9 


Figure  8.  Set-up  for  the  DBD  device  as  described  by  Eq  (47). 

Notice  that  no  free  charge  is  allowed  to  exist  in  the  dielectric,  and  that  the  actual 
numbering  of  the  points  goes  from  0  to  8,  for  a  total  of  9  points,  to  be  consistent  with  the 
formalism  of  starting  array  elements  in  matrices  at  zero.  This  matrix  can  be  easily  solved 
for  by  using  the  Thomas  algorithm  (5:285).  The  Thomas  algorithm  is  a  direct,  single 
pass  method  of  solving  a  tridiagonal  matrix.  The  ID  approximation  used  here  is  justified 
when  the  length  and  width  of  the  electrodes  is  much  greater  than  the  spacing  between 
them,  the  exact  value  of  the  potential  can  then  be  directly  obtained  using  the  Thomas 
algorithm.  This  is  because  the  potential  at  all  points  perpendicular  to  the  point  of  interest 
are  equivalent  in  value  in  this  ID  problem,  which  is  why  the  matrix  shown  in  Eq  (47)  is 
tridiagonal.  However,  in  the  2D  formalism,  an  iterative  scheme  is  used  to  converge  upon 
an  approximately  correct  solution  for  the  scalar  potential  points  over  the  computational 
domain.  This  is  because  the  electrodes  and  free  charges  are  no  longer  considered  to  be 


27 


infinite  sheets  or  over  a  domain  that  is  much,  much  bigger  than  the  spacing  of  the 
electrodes. 

The  iterative  method  invoked  for  the  2D  representation  is  basically  a  5 -point 
averaging  method.  The  four  closest  points  around  the  point  of  interest  are  averaged  to 
find  a  new  value  for  the  fifth  point.  The  notable  difference  between  the  approach  used  in 
this  project  and  the  standard  5 -point  method  is  that  all  the  points  along  a  given  axis  are 
implicitly  solved  for,  using  the  notional  set-up  shown  below.  This  line  method  should 
theoretically  take  a  little  longer  for  each  iteration  of  the  system  than  the  5 -point  averaging 
method,  but  convergence  to  the  set  tolerance  should  happen  in  half  the  number  of 
iterations  (5:285). 


_  Boundary  Conditions 
X  Known  values  at  current  iteration 
.  Values  being  computed 
♦  Values  from  previous  iteration 


Figure  9.  Line  method  used  for  solving  the  2D  potential  profile. 

The  convergence  criterion  of  the  2D  potential  solver  requires  that  all  potential 
points  in  the  computational  domain  vary  by  less  than  1 0  in  value,  between  their 
previous  iteration  value  and  the  current  value.  This  criterion  resulted  in  a  substantial 
number  of  iterations  during  the  first  time  step  but  resulted  in  significant  computational 
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time  savings  for  subsequent  time  step  intervals.  This  is  because  the  line  method 
described  above  saves  and  uses  the  previous  potential  values  from  the  prior  time,  and 
uses  them  as  a  starting  value  for  each  subsequent  time  step.  If  the  time  step  of  the  system 
is  small  enough,  and  the  charged  particles  have  not  moved  very  much,  then  only  a  few 
iterations  are  required  to  find  the  new  potential  values.  For  large  systems,  the  line 
method  described  above  should  give  substantial  computational  time  savings  over  the  5- 
point  method,  but  a  detailed  comparison  was  not  made. 

Values  for  the  electric  field  can  then  be  determined  by  taking  the  negative  of  the 
gradient  of  the  potential  for  every  location  between  points  in  the  computational  area. 
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III.  Computational  Set-up  and  Potential  Solver 


ID  Computational  Set-up 

The  ID  computational  set-up  is  shown  below. 


ID  Dielectric  Barrier  Discharge  Computational  Set-up  for  N=20 


Figure  10.  Characteristic  set-up  geometry  for  the  ID  numeric  model. 

A  staggered  grid  is  employed  where  points  labeled  by  squares  are  where  the 
potential,  number  density,  average  electron  energy  are  solved  for  and  tracked,  as  well  as 
values  for  the  associated  permittivities,  mobilities,  and  diffusion  coefficients.  These 
locations  are  indicated  with  a  whole  number  value.  X’s  indicate  where  electric  field 
values,  particle  and  energy  fluxes,  and  other  half-point  values  identified  as  necessary  for 
the  solution  to  the  three  governing  equations  are  evaluated.  Notice  that  the  dielectric 
boundary  could  extend  over  many  computation  points,  or  simply  be  removed  altogether 
to  simulate  an  exposed  electrode. 

The  first  step  in  determining  whether  or  not  the  numeric  code  used  in  this  project 
was  correct  was  to  validate  the  potential  solver  by  comparing  the  results  to  an  analytically 
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derived  solution  or  to  other  computational  results.  Results  are  shown  below  for  a  ID 


system  with  no  dielectric  barriers  and  no  space  or  surface  charge  present. 


Figure  11.  Potential  profile  for  ID  region  with  no  dielectrics  or  free  charge. 

Here  a  potential  of  100  volts  was  applied  to  the  left  electrode,  while  the  right  electrode 
was  grounded.  The  potential  solver  correctly  returns  a  linear  variation  in  the  absence  of 
any  charge  between  the  electrodes.  Next,  a  dielectric  region  was  added  to  determine  the 
accuracy  of  the  code,  with  various  dielectric  materials  present  between  the  electrodes. 
Shown  are  the  results  for  a  region,  with  a  dielectric  constant  of  4.0,  extending  from  the 
left  electrode  to  midway  between  the  electrodes.  It  is  easily  observed  that  the  value  of 
the  electric  field  (negative  gradient  of  the  potential  profile)  outside  the  dielectric  material 
is  4.0  times  the  value  inside.  This  reduction  in  electric  field  strength  comes  from  the 
polarization  inside  the  dielectric,  caused  by  the  alignment  of  bound  charge  dipoles, 
counteracting  the  applied  field  from  the  electrodes.  This  confirms  that  the  conditions  for 
the  potential  solver  with  dielectric  material  present  are  correct. 
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Figure  12.  Potential  profile  for  inserted  dielectric  and  no  free  charge. 

With  the  addition  of  space  and  surface  charge  to  the  system,  the  computational 
code  for  the  solver  becomes  a  bit  more  complicated,  making  validation  under  these 
circumstances  especially  important.  The  analytic  result  for  the  addition  of  a  uniform 
space  charge  to  a  ID  system,  with  an  electrode  spacing  of  L  is: 
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The  boundary  conditions  at  the  electrodes  are  again  set  to  1 00  volts  and  ground, 
for  the  left  and  right  side  respectively.  In  this  case,  the  uniformly  distributed  positive 
charge  introduced  is  1 0  1 0  C/m  ’ ,  and  the  number  of  evaluation  points  for  the  solver  is 
/V  =  100.  Computational  results,  shown  as  Ws,  are  superimposed  on  the  analytic 
solution,  shown  as  a  line,  in  the  figure  below.  Greater  accuracy  could  have  been  obtained 
at  the  cost  of  computational  time  by  adding  more  points  to  the  simulation;  but  every 
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calculated  point  was  within  10‘6  volts  of  the  analytic  answer,  and  this  was  deemed 
adequate  for  the  project  purposes. 


Figure  13.  Potential  profile  for  no  dielectrics  and  positive  free  charge. 

Finally,  with  the  ID  potential  solver  tested  for  various  dielectric  materials  and 
space  charge,  conditions  needed  to  be  implemented  for  handling  multiple  dielectric 
regions  within  the  computational  area  and  to  account  for  any  accumulated  surface  charge 
on  the  dielectric  surfaces.  The  results  from  these  additions  are  shown  below.  Dielectric 
material  with  a  permittivity  of  4 s0  has  been  introduced  covering  both  electrodes.  Each 
had  a  thickness  of  2  mm,  the  boundaries  of  which  can  been  seen  at  the  X-axis  positions  of 
2  mm  and  8  mm  in  the  figure  below.  Surface  charge  in  the  amount  of  -2.0 -10  6  C/m2 
and  2.0 -10^6  C/m2  has  been  introduced  to  the  left  and  right  dielectric  surfaces 
respectively,  as  well  as  a  uniformly  distributed  positive  charge  between  the  electrodes  of 
1.5-10  4  C/m3 .  It  should  be  noted,  that  Fig  14  below  indicates  that  with  enough  surface 
charge  accumulated  on  the  dielectric  surfaces,  the  electric  field  inside  the  active  region  of 
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a  DBD  device  could  be  significantly  reduced,  leading  to  discharge  extinction  under  DC 
operating  conditions. 

By  defining: 
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where  the  <f>  ’s  represent  potential  values  at  specified  points,  a’s  represent  the  surface 
charge,  p  represents  the  space  charge,  the  s’s  represent  the  dielectric  constant  of  the 
material  and  the  D’s  and  L  represent  the  location  of  the  dielectric  surfaces  and  the 
spacing  of  the  electrodes  respectively,  the  parametric  analytic  solution  of  the  potential  in 
the  three  regions  shown  from  left  to  right  in  Fig  14,  can  be  written: 
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The  location  of  the  pertinent  variables  have  been  labeled  in  the  figure  below.  N  in  the 
calculation  of  the  numeric  solution  has  again  been  set  to  100.  X’s  are  the  numeric 
solution  calculation  points  which  overlay  the  analytic  solution  represented  by  a  line.  The 
two  results  are  in  excellent  agreement,  and  show  that  the  potential  solver  in  ID  is  indeed 
accurate. 


Figure  14.  Potential  profile  for  inserted  dielectrics  with  free  charge. 

With  the  results  of  the  potential  solver  deemed  satisfactory,  the  resultant  potential 
values,  and  thus  electric  field  values  can  be  sequentially  introduced  to  the  continuity  and 
electron  energy  equations  described  in  Chapter  2,  completing  the  last  of  the  unknowns  in 
the  governing  equations.  Shown  below  is  the  highest  level  algorithm  for  describing  the 
iterative  process  by  which  these  equations  will  be  solved  for  in  this  project.  Starting  with 
appropriate  initial  conditions,  the  governing  equations  can  be  solved  in  this  iterative 
manner  to  approximate  the  temporal  evolution  of  the  system.  More  results  of  the  ID 
code  will  be  discussed  in  Chapter  4.  Serious  consideration  was  given  to  using  a 
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completely  implicit  method,  and  solving  for  all  the  equations  simultaneously  at  any  given 
time  step.  However,  the  time  step  size  required  for  stability  reported  by  Marchand 
(15:84)  for  the  fully  implicit  method  was  less  than  that  of  the  sequential  method. 
Therefore,  the  sequential  method  was  chosen  for  both  the  ID  and  2D  geometries.  This 
greatly  simplified  the  numerical  complexity,  while  still  allowing  for  an  accurate  solution 
in  a  reasonable  amount  of  time. 


Figure  15.  Semi-implicit  algorithm  for  solving  the  governing  equations. 
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2D  Computational  Set-up 


The  2D  computational  set-up  is  shown  below. 


Figure  16.  Characteristic  set-up  geometry  for  the  2D  numeric  model. 


Like  the  ID  set-up,  squares  represent  whole  number  values  on  the  grid,  indicating 
where  potential,  number  density  and  electron  energy  density  values  are  tracked  and 
recorded.  X’s  are  at  the  half-point  values,  associated  with  the  fluxes  and  the  electric  field 
values  between  whole  point  locations.  Migrating  the  code  from  a  ID  to  2D  set-up  was 
simplified  by  assuming  a  superposition  principle  for  the  system,  where  the  2D  set-up 
could  be  described  using  its  ID  components.  This  superposition  principle  was 
particularly  advantageous,  given  the  ID  characteristic  of  the  fluxes  in  the  governing 
equations,  as  defined  in  Chapter  2.  Values  for  the  number  and  energy  densities  at  each 
point  were  found  by  solving  the  tridiagonal  governing  equations  along  a  particular  axis. 
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The  changes  in  the  value  of  the  points  due  to  the  flux  from  their  neighboring  points  from 
the  previous  time  step  to  the  current  time  step  were  then  recorded.  This  process  was 
repeated  for  the  other  axis,  using  the  initial  values  from  the  previous  time  step.  The 
changes  due  to  the  X-axis  and  Y -axis  fluxes  were  then  added  to  the  previous  value  to  find 
the  new  value.  This  approach  is  notionally  shown  below. 
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This  piecewise  solving  of  the  new  values  is  incidentally  the  reason  for  the  inclusion  of 
the  Vi  in  the  source  and  loss  terms  of  the  governing  equations  as  discussed  earlier.  This  is 
because  source  and  loss  tenns  are  generally  based  on  the  local  number  density  at  the 
previous  time  step  without  regard  to  direction,  and  are  counted  twice  in  the  superposition 
scheme  implemented. 

Careful  consideration  was  put  into  correctly  describing  the  boundary  conditions 
for  the  2D  environment,  since  the  goal  of  the  project  was  to  model  a  manufacturable, 
aeronautically  realistic  configuration.  The  potential  values  at  the  base  of  the 
computational  area  were  specified  to  be  at  relative  ground,  since  they  would  in  practice 
be  in  contact  with  a  conducting  aircraft  surface.  Two  different  regions  of  variable 
dielectric  constant  were  introduced  to  allow  for  more  flexibility  when  modeling  different 
construction  materials  as  well  as  investigating  the  affects  of  changing  the  insulating 
materials  between  the  aircraft  surface  and  the  electrodes.  The  goal  was  to  correctly 
describe  a  simple  geometry,  which  was  most  like  what  an  experimental  configuration 
might  look  like.  This  project  does  not  implement  code  to  examine  features  of  the  DBD 
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device  if  the  grounded  electrode  is  moved  underneath  the  driven  electrode,  or  to  account 
for  multiple  electrodes  as  in  some  experimental  set-ups  (18:4). 

Charge  buildup  was  allowed  to  accumulate  on  the  top  of  dielectric  #2,  the  affects 
of  which  can  significantly  influence  the  potential  profile  just  as  it  did  in  the  ID  case,  as 
seen  below.  Fig  17  shows  the  potential  profile,  as  well  as  the  resulting  electric  field 
vectors,  for  a  case  in  which  no  surface  charge  is  present.  Fig  18  shows  the  resultant 
potential  and  field  profile  for  a  case  in  which  a  negative  surface  charge  is  unifonnly 
distributed  on  the  dielectric  surface  directly  above  the  grounded  electrode.  The  quantity 
of  the  free  charge  on  the  surface  of  the  dielectric  in  the  second  figure  was  not  recorded, 
and  the  MathCAD  file  used  to  obtain  the  results  was  found  to  be  corrupted.  Therefore, 
the  exact  surface  charge  value  associated  with  Fig  18  is  not  known.  Efforts  to  rebuild  the 
code  were  done,  but  with  no  success.  Completely  rewriting  the  code  was  deemed  too 
time  consuming.  The  presented  results  do  however  serve  as  an  example  of  the  change  in 
electric  field  configuration  under  charge  accumulation  conditions.  The  driven  electrode 
is  held  at  a  potential  of  100  volts.  Notice  the  abrupt  change  in  the  potential  profile  at  the 
dielectric  interface  of  the  Kapton®  and  air.  This  is  because  of  the  factor  of  4.0  change  in 
the  permittivity  of  the  two  mediums.  The  electric  field  vectors  near  the  region  of 
introduced  surface  charge  are  also  significantly  different  between  the  two  figures, 
suggesting  that  the  accumulated  free  charge  on  the  surface  of  the  dielectric  can  play  an 
important  role  in  the  dynamics  of  the  system,  and  therefore  the  flow  control  associated 
with  various  operating  parameters. 
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X-axis  (arbitrary  units) 


Figure  17.  2D  potential  profile  and  associated  E  field  vectors  for  no  free  charge. 
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X-axis  (arbitrary  units) 

Figure  18.  2D  potential  profile  and  associated  E  field  vectors  for  negative  surface  charge. 
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One  more  potential  solver  verification  was  made  for  the  2D  code  before 
integrating  it  with  the  governing  equations  to  solve  the  dynamics  of  the  system.  The 
value  of  the  normal  component  of  the  displacement  ( D± )  at  the  interface  of  the  dielectric 
and  gas  of  the  system  should  be  continuous  and  approach  the  analytic  solution  of: 

^Ldieletric  ^ Xfreespace  *  (^^0 

It  should  be  noted  that  the  individual  X-axis  and  Y-axis  component  values  of  the 
displacement  can  be  known  for  the  aforementioned  2D  geometry  because  of  the  systems 
boundary  conditions,  and  the  iterative  method  used  to  determine  the  potential  values. 
This  allows  for  the  component  contribution  to  the  divergence  of  the  displacement  by  any 
free  charge  to  be  found.  If  the  specified  boundaries  and  spatial  discretization  of  the 
system  are  adequate,  then  one  would  expect  that  the  perpendicular  component  of  the 
displacement  vector  of  the  solved  system  approach  or  come  close  to  the  real  value. 

For  this  case,  no  surface  charge  will  be  allowed.  Rewriting  the  Eq  (50)  above  into 
its  equivalent  dielectric  constant  and  potential  components  representation  along  the  Y- 
axis  gives: 


'  dieletric 


(tx.y  ~  1  ) 

Ax 


=  e 


freespace 


(0x.y  + 1  ~  A.y) 

Ax 


(51) 


since  the  values  of  the  displacement  vector  are  solved  for  at  locations  of  ±  —  Ax  above 

and  below  the  interface  point.  By  further  manipulating  the  equation  above,  it  follows 
that: 
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dieletric 


(52) 


£ 


£ 


freespace 


Ax _ 

(</>,, y  ~  #r,y- 1  ) 

Ax 


If  a  nominal  values  of  £f.eespace  ~  1-0  and  sdielectric  ~  4.0  are  specified  for  the  potential  solver, 


it  follows  that  both  sides  of  the  equation  above  should  evaluate  to  4.0  in  the  limit  that  the 
grid  points  become  arbitrarily  close  to  the  dielectric  interface.  The  convergence  to  the 
value  4.0,  for  a  point  on  the  interface,  lying  halfway  between  the  start  and  end  of  the 
grounded  electrode  for  increasing  values  of  N,  corresponding  to  smaller  and  smaller 
spacing  between  the  potential  points,  is  shown  in  the  table  below.  The  length  and  width 
of  the  computational  region  is  0.1  meter,  defining  Ax  =  0.1/iV,  and  the  convergence 
tolerance  was  set  tolO  lS.  Values  of  the  time  to  convergence  for  the  solver  were  also 
recorded,  as  well  as  their  associated  computational  iteration  cost.  Notice  that  for 
increasing  value  of  N  (or  more  closely  spaced  calculated  potential  points,  given  better 
spatial  resolution  to  the  system)  the  time  to  convergence  is  approximately  quadrupled  for 
every  doubling  of  the  resolution.  The  potential  solution  is  in  fact  the  most 
computationally  intensive  portion  of  the  total  solution,  so  the  number  of  calculation 
points  was  held  at  N  =  50  for  the  remainder  of  the  project,  to  allow  for  results  to  be 
obtained  within  a  reasonable  amount  of  time. 


Table  1.  Computational  values  for  the  electric  field  converging  on  the  analytic  solution. 


Value  of  N 

F  F 

L Above  j  1. Below 

Iterations  to  Convergence 

Time  to  Convergence 

50 

3.312 

~  8,000 

<30  sec 

100 

3.742 

~  32,000 

~  1  min 

150 

3.840 

~  74,000 

~  3  min 

200 

3.862 

*  124,000 

~  10  min 

250 

3.883 

~  190,000 

~  45  min 
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The  2D  potential  solver  was  further  validated  by  coding  its  geometry  to  simulate 
the  ID  environment.  This  was  done  by  modifying  the  boundary  conditions  to  make  the 
top  and  bottom  of  the  computational  region  electrodes.  The  side  boundary  conditions 
remained  set  to  ensure  no  perpendicular  field  components,  thus  reproducing  the  ID 
approximations.  The  results  from  the  modified  2D  code,  matched  the  ID  code  perfectly 
for  a  series  of  test,  including  single  dielectric,  multiple  dielectrics,  and  space  and  surface 
charge  presence,  and  the  analytically  derived  uniform  space  charge  presence.  A 
representative  result  of  this  pseudo  ID  code  is  shown  in  Fig  19  below.  No  free  charge 
was  present  in  the  system  and  a  dielectric  of  relative  constant  4.0  was  inserted  from  the 
left  electrode  to  midway  between  the  electrodes.  Comparing  this  2D  result  with  the  ID 
result  shown  in  Fig  11,  shows  that  they  are  identical.  The  last  validation  for  the  2D  code 
was  done  by  using  the  standard  5 -point  averaging  method  in  Microsoft  Excel  under 
identical  point  geometry  to  solve  for  the  potential  profile  without  space  charge  present, 
and  comparing  the  outputs  of  both  programs.  As  seen  below  in  Fig  20  for  a  small  portion 
of  the  results  near  the  electrodes,  the  comparison  matched  up  perfectly,  and  helped 
established  further  confidence  in  the  2D  potential  solver. 


2D  Numeric  Result  Simulating  (tie  ID  Environment  for  Various  Perspectives 


Figure  19.  The  2D  potential  solver  results  for  a  pseudo  ID  geometry. 
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Excel  Results 


41 .5892644 

38.64189843 

35.35925385 

44.543817119754 

40.273882505125 

35.883112142282 

50 

42.026702330137 

36.004220219998 

50 

41 .828706595423 

34.975244717763 

50 

40.312879333792 

32.583516441306 

50 

3B  346187571 105 

36.839294298439 

28.685924142656 

28.358373717306 
21 .556090243976 

26.797795318846 

20.002124457105 

12.645952478869 

19.422864582977 

1 1 .878825888048 

0 

14.474144191556 

8.0903145121 1 1 

0 

11.129810128911 

6.008287968841 

0 

8.804726706153 

4.813027234343 

0 

7.149608047216 

4.439094262376 

1 .912694719697 

5.711456839952 

3.881047048250 

2.283872976059 

4. 158926e+001 
4 . 454B82e+001 
5 . OOOOOOe+OOl 
5 . OOOOOOe+OOl 
5. OOOOOOe+OOl 
5. OOOOOOe+OOl 
3 . 634619e+001 
2 . 679780e+QQl 
1 . 942286e+001 
1 . 447414e+001 
1 . 112981e+001 
8 . 804727e+000 
7 . 149608e+000 
5 . 711457e+000 


Numeric  Code  Results 
3 . 864190e+001 
4 . 027388e+001 
4 . 202670e+001 
4. 182871e+001 
4. 031288e+001 
3 . 683929e+001 
2 . 868592e+001 
2 . 000212e+001 
1 . 187883e+001 
8. 090315e+000 
6. 008288e+000 
4. 813027e+000 
4.439094e+000 
3 . 881Q47e+G0Q 


3 . 5  3  592  5  e+OOl 
3. 588311e+001 
3 . 600422e+001 
3 . 497524e+001 

3 . 2  583  52e+001 
2 . 83  5837e+001 
2 . 15  5609e+001 

1. 2  64  595  e+OOl 
0 . OOOOOOe+OOO 
0 . OOOOOOe+OOO 
0. OOOOOOe+OOO 
0. OOOOOOe+OOO 
1.912695e+000 

2 . 2  8  3  87  3  e+000 


Figure  20.  5-point  averaging  Excel  and  numeric  code  results  for  similar  regions. 


44 


IV.  Results  and  Conclusions 

The  results  for  various  cases  using  the  ID  code  are  shown  and  discussed  below. 
Electrode  potentials  are  modeled  as  DC  for  all  the  results,  even  though  DBD  actuators  are 
run  in  an  AC  mode.  The  DC  approximation  is  sufficient  for  the  time  scales  investigated. 
The  time  step  used  for  all  the  simulations  was  10  10  seconds.  This  is  about  two  orders  of 
magnitude  smaller  than  what  the  dielectric  relaxation  time  tdi ,  defined  by: 

tdi  =  7  I  f°  ,  (53) 

q(n  /j  +n  /j  ) 

was  expected  to  be,  but  still  not  so  small  as  to  require  an  inordinate  amount  of 
computational  time  per  simulation. 

The  code  was  migrated  to  the  2D  environment;  but  unfortunately  inherent 
problems  with  the  boundary  conditions  in  the  governing  equations,  instabilities  and  errors 
in  the  code,  discovered  late  in  the  project,  kept  the  ultimate  goal  of  modeling  a  realistic 
interaction  between  the  charged  species  and  the  neutral  particles  from  being 
accomplished. 

DC  Discharge  Sheath  Test 

The  first  characteristic  test  chosen  to  see  if  the  electron  and  ion  continuity 
equations  returned  accurate,  expected  results  when  integrated  with  the  potential  solver 
and  solved  sequentially,  was  to  look  at  the  DC  discharge  sheath  that  forms  during  the  first 
few  nanoseconds  of  the  discharge.  Under  DC  operating  conditions,  it  should  be  expected 
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that  the  electrons,  because  of  their  relatively  high  mobility,  as  compared  to  the  ions 
would  move  under  any  applied  field  in  the  system  quickly,  leaving  the  resulting  plasma 
with  a  sheath  formation.  This  sheath  of  positively  charged  plasma  should  remain  close  to 
the  cathode  as  electrons  are  forced  toward  the  anode.  This  depletion  causes  a  secondary 
electric  field,  that  acts  to  shield  the  rest  of  the  plasma  from  the  anode,  and  it  remains  in  its 
quasi-neutral  state.  A  sheath  region  should  also  develop  around  the  anode,  as  electrons 
are  lost  to  the  surface  of  the  electrode  or  dielectric  coating  due  to  the  applied  electric  field 
and  their  high  thermal  velocities.  Results  for  the  ID  system  with  (Fig’s  21,  22,  and  23.) 
and  without  (Fig’s  24,  25,  and  26.)  dielectrics  as  thin  coatings  on  the  electrode  surfaces 
are  shown  below.  Charge  densities  for  the  electrons  and  ions  are  plotted  normalized  to 
their  initial  condition,  which  in  these  cases  was  a  uniform  distribution.  Parameters  for 
this  test,  with  Ar  gas  and  the  anode  held  constant  at  -50  volts,  were  held  constant  as 
follows: 

ne  =  n‘  =  101!  /m3  u‘ =  1.0  eF  P  =  \00Vorr 

/j‘  =0.3  m2  IV  Is  ju‘  =  \0^m2  IV  Is  v  =  0.0 

D‘  =0.3  m2  Is  D‘  =  I0~4m2  Is  0  =  0.0 
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Distance  between  Electrodes  (in  meters) 


Figure  2 1 .  Potential  profile  during  formation  of  the  DC  discharge  sheath  in  a  system  with 
bare  electrodes. 


Distance  between  Electrodes  (in  meters) 


Figure  22.  Normalized  electron  density  during  formation  of  the  DC  discharge  sheath  in  a 
system  with  bare  electrodes. 


Distance  between  Electrodes  (in  meters) 


Figure  23.  Normalized  ion  density  during  formation  of  the  DC  discharge  sheath  in  a 
system  with  bare  electrodes. 
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The  width  of  the  sheath  region  for  the  results  above  is  approximately  0.001  m. 
This  is  consistent  with  the  method  used  by  Poggie  to  estimate  the  region  width  (16:2). 
The  characteristic  width  should  be  on  the  order  of  a  Debye  length  defined  as  follows: 

1 

^ Debye 

Estimating  the  electron  temperature  to  be  approximately  10,000°  K  or  about  1.0  eV,  it 
can  be  shown  that  ADebve  =  0.000218  m ,  indicating  that  the  sheath  region  observed  in 

these  calculations  is  about  five  Debye  lengths  wide.  This  is  of  the  order  that  was 
expected,  and  was  therefore  deemed  a  reasonable  result. 

Notice  that  the  charge  accumulated  on  the  dielectric  surfaces  of  the  figures  below 
affects  the  potential  profile  and  width  of  the  sheath  region.  It  is  also  seen  that  many  more 
charged  particles  are  left  in  the  system  with  dielectric  coatings.  This  is  because  of  the 
reduced  field  caused  by  the  accumulation  of  charge  on  the  dielectrics,  and  the  fact  that 
the  charges  are  not  allowed  to  leave  the  system  through  the  electrode  as  before. 


e  2 

n  a 


(54) 
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Distance  between  Electrodes  (in  meters) 


Figure  24.  Potential  profile  during  formation  of  the  DC  discharge  sheath  in  a  system  with 
dielectrics  coatings. 


Distance  between  Electrodes  (in  meters) 


Figure  25.  Normalized  electron  density  during  formation  of  the  DC  discharge  sheath  in  a 
system  with  dielectric  coatings. 


Distance  between  Electrodes  (in  meters) 


Figure  26.  Normalized  ion  density  during  formation  of  the  DC  discharge  sheath  in  a 
system  with  dielectric  coatings. 
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The  results  showed  good  agreement  of  the  electron  density  with  previous  numeric 
calculations  as  seen  in  Fig  27.  In  Fig  27  the  number  density  is  not  normalized  to  the 
initial  condition.  The  ion  densities  differ,  and  this  subsequently  affects  the  calculated 
potential  profile  (12:11). 


Figure  27.  DC  discharge  results  reported  by  Hilbun  (12: 1 1) 

Ambipolar  Diffusion 

For  no  applied  potential  across  the  electrodes  and  an  initially  quasi-neutral  plasma 
distribution,  the  charged  particle  density  near  the  electrodes  should  exhibit  a  rapid 
depletion  of  the  electrons,  due  to  their  much  higher  thermal  speed  than  that  of  the  ions. 
As  electrons  are  lost  to  the  boundaries  much  more  quickly  than  the  ions,  an  induced 
electric  field  is  formed  that  tries  to  bring  the  system  back  to  its  quasi-neutral  state  by 
balancing  the  fluxes.  This  process  of  the  plasma,  where  diffusion  dominates  the  first  few 
time  steps  and  then  the  induced  field  effect  starts  to  limit  the  process,  is  seen  in  Fig’s  28 
and  29  for  various  time  steps  below.  The  initial  conditions  for  the  system  were  that  both 
electrodes  be  held  at  ground  with  dielectric  coatings  present,  and  that  the  plasma  be 
distributed  sinusoidally  as  follows: 
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ne  =  nx  =  n(0)Sin 


^  nx^ 

L 


(55) 


with  /?(())  =  1014  /m3.  The  transport  coefficients  were  determined  using  the  BOLSIG  fits 
with  the  electron  energy  held  constant  at  1.0  eV,  in  a  1  Torr  gas  of  N2.  Ionization  and 
recombination  were  set  to  0.  The  plots  below  in  Fig’s  28  and  29  show  the  relative 
changes  in  number  density  as  the  system  transitions  from  a  diffusion  to  an  ambipolar 
particle  flux. 
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Figure  28.  Electron  density  during  ambipolar  diffusion. 
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Figure  29.  Ion  density  during  ambipolar  diffusion. 
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This  non-uniform  decay  shown  above  is  indicative  of  the  induced  electric  field  which 
results  from  the  initial  charge  separation  causes  by  the  fast  moving  electrons.  It  should 
also  be  noted  that  the  system  has  not  yet  reached  an  equilibrium  point  for  the  time  steps 
shown.  Fig  29  shows  an  initial  build-up  of  positively  charged  ions  near  the  dielectric 
boundaries,  and  could  be  indicative  of  a  bottleneck  in  the  numeric  method  due  to  an 
unrealistic  boundary  conditions,  or  an  error  in  the  code.  The  source  of  the  apparent 
problem  was  never  clarified.  An  immediate  rise  in  the  potential  profile  over  the  entire 
domain  was  also  noticed  as  seen  below;  but  then  began  to  settle  back  to  its  initial 
condition,  as  should  be  expected  as  the  system  tries  to  return  to  a  quasi-neutral  state. 
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Figure  30.  Potential  profile  during  ambipolar  diffusion. 


Furthermore,  the  accumulated  flux,  calculated  by  the  numeric  code  at  the 
dielectric  surfaces  (which  are  in  this  case  one  Ax  thick,  and  placed  directly  over  both  the 
anode  and  the  cathode),  should  average  to  zero  after  the  system  reaches  a  steady  state 
distribution,  as  well  as  be  nominally  equal  to  each  other  as  the  dynamics  of  the  system 
develop  because  of  the  symmetry  involved.  Shown  in  the  table  below  are  the  results  for 
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the  surface  charge  per  meter  squared  at  both  dielectric  surfaces  for  various  times.  Note 


that  the  values  correctly  approach  zero  as  time  increases. 


Table  2.  List  of  calculated  accumulation  of  surface  charge  on  the  dielectric  surfaces  for 
various  times  during  ambipolar  diffusion. 


Time  =  0  sec 

Time  =  lgs 

Time  =  2gs 

Time  =  3gs 

Time  =  4gs 

Time  =  1ms 

Surface  Charge  on  left 
dielectric 

0 

-2.0421e-7 

-1.5036e-7 

-5.3790e-8 

-3.9940e-8 

-2.2720e-10 

Surface  Charge  on  right 
dielectric 

0 

-2.0406e-7 

-1.5025e-7 

-5.3680e-8 

-3.9830e-8 

-1.2580e-10 

Constant  Characteristic  Parameters  Comparison 

With  the  continuity  equations  tested,  a  full  implementation  with  the  electron 
energy  density  equation  included,  as  well  as  the  associated  rate  and  transport  fit,  was 
tested.  This  final  test  was  conducted  on  the  ID  code  on  a  system  comprised  of  100%  Nt 
by  changing  the  gas  pressure,  electrode  spacing,  and  initial  charge  densities  for  systems 
that  held  the  parameters  below  constant,  and  observing  their  behavior  over  time: 


N 


—  ,  Pxd,  n 


5,0 


Neutral 


N 


Neutral 


The  first  tenn  is  called  the  reduced  field,  the  second  is  the  pressure  multiplied  by  the 
electrode  spacing,  commonly  associated  with  the  Paschen  curve  (which  describes  the 
breakdown  voltage  relationship  of  the  system  as  a  function  of  the  pressure  of  the  system 
and  the  electrode  spacing),  and  the  final  term  is  the  fractional  ionization  of  the  system. 
E  is  the  electric  field,  ns'°  is  the  initial  uniform  charge  density  for  both  species,  d  is  the 
electrode  separation  distance,  and  P  the  pressure  of  the  neutral  gas  which  is  related  to 
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N Neutml-  This  neutral  number  density  was  previously  defined  in  Eq  (15).  The 
temperature  of  the  system  is  set  to  300°  K.  The  driven  electrode,  the  anode  in  this  test 
case  is  located  on  the  left  of  the  ID  environment.  The  potential  of  the  anode  was  fixed  at 
2.5  kV  for  all  cases.  This  occurs  naturally  by  holding  the  reduced  field  term  defined 
above  as  constant.  The  governing  equations  retain  the  same  form  for  all  cases,  but  values 
for  the  neutral  number  density,  the  electric  field,  and  the  rate  and  transport  coefficients 
change.  The  initial  rate  and  transport  coefficients  remained  the  same  however,  because 
the  initial  value  given  to  the  local  electron  energy  was  defined  as  3.0  eV  at  every  location 
in  the  system.  These  changes  resulted  in  different  behaviors  being  incorrectly  observed 
for  each  simulated  scenario  in  this  project,  in  part  because  of  the  limitations  associated 
with  the  rate  and  transport  coefficients  as  well  as  the  spatial  discretization,  which 
changed  as  the  electrode  separation  changed.  In  reality,  the  results  should  have  been 
identical. 


Table  3.  Values  used  for  comparison  between  constant  characteristic  systems. 


Pressure  (in  Torr) 

Electrode  Separation  (in  m) 

Charge  Density  (in  m-3) 

Potential  (in  V) 

P=1 

d=.l 

ns’°=10ls 

V=2500 

P=10 

d=.01 

ns-°=1016 

V=2500 

'"d 

II 

o 

o 

d=.001 

ns-°=1017 

V=2500 

P=350 

d=2.857xl0'4 

ns'°=3.5xl()17 

V=2500 

P=760 

d=1.316xl0'4 

ns'°=7.6x  1 01 7 

V=2500 

The  goal  of  this  test  was  to  compare  these  different  conditions  for  various  times, 
and  report  the  resultant  charge  and  energy  densities;  as  well  as  the  relative  dynamics  of 
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the  system.  During  a  discussion  with  Bailey  (1)  it  was  discovered  that  using  a 
transformation  of  parameters  approach,  it  was  possible  to  fully  describe  the  results  of  the 
simulation,  for  any  combination  of  the  characteristic  parameters  listed  above,  given  a 
solution  under  any  other  conditions.  This  would  have  been  advantageous  when 
comparing  any  previously  published  results,  as  well  as  analytic  solutions  if  available. 
However,  as  seen  in  the  figures  below,  unstable  behavior  was  observed.  The  charge 
densities  of  the  P  =  1 00  Torr  and  greater  cases  completely  disappeared  during  the  first 
microsecond  of  operation.  For  simulation  run  times  on  the  order  of  nanoseconds, 
unpredictable  instabilities  at  these  higher  pressures  were  observed,  but  not  reported  on 
here.  These  instabilities  were  thought  to  arise  from  the  values  returned  by  the  electron 
energy  equation  solution.  These  energy  density  values  exhibited  extreme  oscillatory 
behavior,  and  negative  values  were  sometimes  generated.  The  location  of  the  problems 
appeared  to  occur  at  the  boundaries  of  the  system,  but  was  found  to  occur  at  other 
locations  as  well  for  longer  run  times.  Therefore,  the  boundary  conditions  of  the  system 
were  determined  to  not  be  the  cause  of  the  instability.  A  good  deal  of  work  went  into 
isolating  the  problems  with  the  energy  equation;  but  it  was  eventually  eliminated  from 
the  program,  and  replaced  by  another  method  for  determining  the  local  rate  and  transport 
values  as  discussed  in  the  next  section. 

In  Fig  3 1  below,  it  is  seen  that  the  electron  density  behaves  as  discussed  in  the 
DC  discharge  sheath  section  for  the  10  Torr  case  near  the  electrodes;  but  exhibits 
instabilities  which  eventually  return  to  a  steady  state  value,  in  the  central  portion  of  the 
discharge  as  seen  in  Fig  32  and  Fig  33.  An  initial  build-up  of  electrons  near  the  anode  in 
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the  1  Torr  case  in  Fig  31,  not  present  at  later  times,  is  seen  and  gives  suspicion  of  an 
underlying  numeric  error.  The  results  in  Fig  31  look  to  be  very  similar  to  a  steady  state 
solution  during  the  first  microsecond  of  operation,  and  should  have  been  investigated  for 
earlier  times.  A  magnified  view  (Fig  33)  over  various  times  for  the  instabilities  shown  in 
Fig  32  is  presented. 


Normalized  Distance  Between  Electrodes 

-  P=1  Torr 

-  p=10  Torr 


Figure  3 1 .  Charge  distribution  after  1  /us  for  constant  characteristics  comparison. 


Normalized  Distance  Between  Electrodes 

-  P=1  Torr 

-  P=10  Torr 


Figure  32.  Charge  distribution  after  3  /us  for  constant  characteristics  comparison. 
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Figure  33.  Magnified  view  of  the  instabilities  near  the  anode  encountered  at  P=10  Torr. 


ID  Fixes 


Eventually,  the  erratic  behavior  of  the  numeric  code  was  linked  to  the  values  of 
the  local  average  electron  energy  which  are  used  to  determine  the  rate  and  transport 
coefficients  as  previously  described.  The  setting  of  minimum  and  maximum  values  of 
these  energies  in  an  attempt  to  keep  the  results  within  the  range  of  goodness  for  the  fits, 
masked  a  problem  with  the  solution  of  the  energy  equation.  Negative  numbers,  which 
began  appearing  at  the  boundaries  and  sometimes  at  apparently  random  locations 
throughout  the  computational  domain,  created  serious  problems  for  the  numeric  code. 
Therefore  the  energy  equation  was  eliminated,  and  replaced  by  the  reduced  field 
representation  for  determining  the  rate  and  transport  coefficients.  This  reduced  field 
representation,  previously  defined  as: 

E 

Reduced  Field  = - , 
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is  a  characteristic  of  the  conditions  in  the  plasma,  and  is  related  to  the  local  average 
electron  energy.  Instead  of  having  to  rewrite  the  rate  and  transport  coefficient  fit  portion 
of  the  numeric  code,  the  BOLSIG  solver  was  once  again  used  to  relate  the  reduced  field 
to  the  average  local  electron  energy  as  seen  below.  A  fit  of  this  relationship  was 
determined,  and  used  in  conjunction  with  the  previous  fits  to  assign  new  coefficients  to 
all  points  in  the  system  as  seen  in  Fig  34  below.  This  approach  also  had  difficulties 
associated  with  it,  and  the  range  over  which  the  coefficients  could  be  accurately 
described  was  further  reduced  because  of  the  goodness  of  the  fits.  This  limited  the  range 
of  input  parameters,  such  as  electrode  potential  values  and  pressures  to  a  specific  range. 
Unfortunately,  this  range  did  not  include  the  specified  parameter  values  for  the  constant 
characteristics  comparison  test  discussed  above.  The  initial  guess  of  3.0  eV  for  the 
average  electron  energy  for  this  test  case  reported  above  was  found  to  be  in  error,  and 
was  most  likely  the  main  contributor  to  the  erratic  behavior  observed.  This  value  did  not 
allow  for  accurate  initial  rate  and  transport  coefficients  and  limited  the  ability  of  the 
simulation  to  produce  self  sustaining  ionization.  The  correct  value  for  the  initial  average 
local  electron  energy  value  should  have  been  approximately  13  eV,  as  estimated  using  the 
reduced  field  characteristic.  This  energy  value  is  well  outside  the  range  for  which  the 
rate  and  transport  fits  were  calculated,  and  presents  a  significant  limitation  to  the 
suitability  of  this  numeric  code  for  modeling  realistic  systems.  Unfortunately,  this 
oversight  was  not  found  until  very  late  in  the  project,  and  time  constraints  did  not  allow 
for  it  to  be  corrected.  New  BOLSIG  fits  over  a  wider  range  of  electron  energies  or  direct 
reduced  field  value  fits  introduced  into  the  code  may  take  care  of  this  limitation  in  the 
future. 
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N2  BOLSIGData-  Reduced  Field  vs.  Energy 


Figure  34.  BOLSIG  fits  for  relating  the  reduced  field  to  the  local  average  electron  energy. 

With  the  electron  energy  equation  eliminated  and  the  conditions  under  which  the 
code  could  be  run  more  definitively  specified,  it  was  found  that  the  source  of  the  negative 
energy  values  was  actually  from  the  solution  of  the  electron  continuity  equation.  Under 
certain  time  step  and  computational  cell  spacing  conditions,  it  was  found  that  the  solution 
to  the  governing  equations,  especially  the  electron  continuity  equation,  gave  inaccurate 
results  at  the  boundaries  because  of  the  interaction  of  the  thermal  velocity  used  to 
describe  the  boundary  flux,  and  the  local  field  and  mobility  coefficient.  Special  care  had 
to  be  taken  to  not  let  the  rate  at  which  the  electrons  moved  through  the  system,  exceed 
what  the  spatial  discretization  of  the  model  could  accurately  handle.  This  was  done  by 
carefully  choosing  the  distance  between  cells  and  the  time  step.  By  using  the  expected 
thermal  velocities  of  the  electrons  under  the  conditions  tested,  the  cell  spacing  was 
adjusted  so  that  the  electrons  were  not  allowed  to  move  more  than  %  of  a  cell  during  each 
time  step.  As  seen  in  Eq  (56)  below  (which  is  Eq  (18)  rewritten  for  clarity)  the  method 
used  to  solve  the  governing  equations,  could  potentially  overstep  its  bounds,  leading  to  an 
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erroneous  negative  result  for  the  new  time.  This  occurred  when  the  field  driven  flux  and 


the  thermal  flux  were  oppositely  directed. 
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Various  instabilities  were  still  found  after  all  these  fixes  had  been  implemented, 
and  no  correlation  between  when  or  where  they  would  present  themselves  was  found. 
This  made  it  difficult  to  focus  in  on  any  other  areas  to  fix  problems.  No  further  ID 
results  where  obtained. 


2D  Results  and  Problems 

The  ID  and  2D  codes  were  developed  simultaneously,  and  therefore  problems 
found  in  one  had  to  be  addressed  in  the  other.  Obtaining  even  preliminary  results  from 
the  2D  code  proved  difficult  because  of  a  lack  of  suitable  initial  conditions.  The  plan  was 
to  introduce  a  Gaussian  distribution  of  charge  as  seen  in  Fig  35  below,  disallow  any 
charge  buildup  on  the  dielectric  surface,  and  run  the  simulation  until  a  steady  state 
solution  was  reached.  The  exposed  driven  electrode  can  be  seen  towards  the  bottom  left 
of  the  domain,  while  the  grounded  electrode  is  buried  in  the  dielectric  material.  The  peak 
charge  density  seen  below  is  about  1015  m  3 ,  which  falls  off  rapidly  to  zero  in  the  outer 
regions  of  the  computational  area. 
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Figure  35.  Gaussian  used  as  an  initial  condition  for  charge  density  distribution. 

Due  to  the  previously  described  difficulties  and  problems  finding  erratically 
behaving  cells  in  the  much  larger  2D  geometry,  the  convergence  criteria  was  never  met, 
and  therefore  no  suitable  initial  conditions  could  be  implemented.  Definite  changes  in 
the  density  profile,  charge  accumulation,  and  evolution  of  the  potential  profile  over  time 
were  observed  and  looked  nominally  like  one  would  expect,  as  seen  below.  Here  the 
dimensions  of  the  2D  system  reflect  a  physically  realistic  set-up,  with  the  base  of  the 
computational  area  grounded,  the  driven  electrode  in  the  bottom  left  of  the  region  is  held 
at  2.5kV,  and  the  grounded  electrode  is  towards  the  bottom  right.  Rate  and  transport 
coefficients  are  determined  using  the  reduced  field  representation  fix  from  the  ID 
environment,  and  the  initial  density  distribution  is  the  Gaussian  from  Fig  35  above.  The 
100%  iV2  background  gas  pressure  is  1  Torr,  and  the  driven  electrode  is  exposed. 

Secondary  emission  of  electrons  due  to  electrode  impact  by  ions  is  turned  on,  but  not 
expected  for  this  configuration  since  the  driven  electrode  is  held  at  a  positive  potential. 
The  contour  lines  for  the  resultant  potential  profile  over  time  are  shown  in  Fig  36,  as  well 
as  the  evolution  from  the  initial  Gaussian  distribution  of  the  electron  number  density  in 
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Fig  37.  Notice  the  erratic  behavior  of  the  density  on  the  left  side  of  the  computational 
region  around  the  driven  electrode  in  Fig  37.  A  suitable  explanation  for  this  behavior 
was  not  found.  Also  notice  the  charge  buildup  on  the  dielectric  surface  in  Fig  38.  It 
progressed  from  an  initial  condition  of  zero  everywhere,  and  shows  the  greatest  buildup 
in  regions  closest  to  the  electrodes.  Overlaid  is  the  representative  initial  electron  density 
distribution  at  the  points  directly  above  the  dielectric  surface.  The  shape  of  which  seems 
to  correspond  to  the  calculated  charge  accumulation,  but  notice  the  erratic  behavior  of  the 
charge  density  to  the  left  of  the  electrode.  The  value  of  the  reported  charge  accumulation 
points  can  be  estimated  by  considering  the  initial  conditions  of  the  electron  thermal 
velocity,  the  boundary  conditions,  and  the  electron  number  density.  This  estimate 
showed  good  agreement  to  the  calculated  value,  so  the  charge  accumulation  was 
considered  in  line  with  what  should  be  expected.  The  charge  buildup  is  thought  to 
migrate  down  the  length  of  the  buried  electrode  as  time  progresses,  and  can  narrowly  be 
seen  upon  close  inspection  and  interpretation  of  the  potential  profile  plots  in  Fig  36 
below. 


No  confidence  was  placed  in  these  results  however  because  of  the  difficulties  in 
reaching  a  steady  state  solution  associated  with  the  2D  code  that  were  encountered. 
These  problems  were  indicative  of  either  program  coding  errors  or  basic  physical  model 
assumptions  and  techniques  that  need  to  be  worked  out  before  any  more  work  is  done. 
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Figure  36.  Plots  of  the  2D  potential  profile  progression  for  various  times.  Changes  in  the 
profile  are  seen  as  charge  moves  in  the  system  and  accumulates  on  the  dielectric  surface. 
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Electron  Density  at  t=40  ns 
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Figure  37.  Electron  number  density  for  various  times  after  initial  Gaussian.  Instabilities 
around  the  exposed  electrode  area  are  seen. 
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Charge  Accumulation  oil  Dielectric  Surface  after  160  ns 
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Figure  38.  Charge  buildup  on  the  dielectric  surface  at  end  of  simulation.  The  region 
between  7  and  17  units  on  the  axis  is  the  location  of  the  exposed  electrode.  Therefore,  no 
charge  was  allowed  to  accumulate  there. 
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V.  Summary  of  Work  and  Recommendations  for  Future  Efforts 
Summary  of  Work 

This  project  developed  three  governing  equations  used  to  characterize  the 
macroscopic  properties  of  the  plasma  dynamics  of  a  ID  and  2D  plasma  actuator  set-up. 
Methods  for  determining  the  rate  and  transport  coefficients,  solving  for  the  potential 
profile  of  the  system,  and  integrating  all  the  pieces  into  a  sequential  solving  algorithm 
were  presented.  Preliminary  results  were  found  for  the  DC  discharge  sheath  model,  the 
ambipolar  diffusion  case,  and  by  holding  some  characteristic  parameters  of  the  system 
constant.  The  main  goal  of  developing  an  accurate  predictor  tool,  in  good  agreement  to 
the  experimentally  observed  flow  control  and  boundary  layer  reattachment  forces,  for 
determining  the  momentum  transfer  to  the  neutral  air  flow  of  a  system  by  the  moving 
plasma,  was  not  accomplished.  However  an  accurate  potential  solver  in  both  ID  and  2D 
was  shown  and  that  represents  a  suitable  achievement. 

Recommendations  for  Future  Efforts 

Various  improvements  can  be  implemented  to  the  current  code  in  the  future. 
Changes  to  the  dimensions  of  the  computational  area  should  be  made  in  order  to  save 
computational  time  in  regions  of  low  activity,  as  well  as  allowing  for  the  simulation  of 
more  realistic  rectangular  like  regions,  with  many  electrodes  present  (as  seen  from 
experimental  results  in  Fig  39  below),  and  a  higher  resolution  view  of  the  progression  of 
the  plasma  close  to  the  boundary  surface.  The  placement  of  grounded  and  driven 
electrodes  should  be  examined,  and  allowed  to  overlap  to  study  if  any  benefit  can  be 
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obtained  through  this  geometry.  Grid  spacing  should  be  allowed  to  vary  from  place  to 
place  within  the  domain,  to  more  accurately  model  the  interfaces  and  the  geometry  of  the 
electrodes  in  the  2D  case.  Better  characterization  of  the  dielectric  surface  where  charge 
build-up  is  allowed  to  either  flow  through  the  dielectric  material  due  to  the  applied  field 
or  jump-off  to  contribute  to  the  ionization  under  reversing  field  conditions  should  be 
looked  at.  The  boundary  conditions,  and  the  cause  of  the  instabilities  present  in  the  ID 
and  2D  simulations  described  above  also  need  more  examination. 


The  most  computationally  intensive  portion  of  the  program  is  the  potential  solver, 
and  it  would  be  very  interesting  to  look  at  other  methods  for  solving  the  potential, 
especially  in  2D.  By  looking  at  the  traditional  5-point  method  it  can  be  seen  that  the  full 
matrix  for  solving  the  system  implicitly  has  a  pentadiagonal  like  symmetry  as  seen  below 
in  Fig  40.  The  uppermost  and  lowermost  diagonals  vary  in  their  position  depending  on 
the  size  of  the  array,  but  the  center  tridiagonal  feature  remains.  The  figure  below  shows 
the  basic  set-up  for  a  6X6  grid  where  the  bottom  row  is  grounded,  and  the  sides  and  top 
are  limited  by  boundary  conditions  to  make  the  perpendicular  component  of  the  electric 
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field  equal  to  zero.  As  before,  this  condition  is  physically  unrealistic,  except  in  cases 
where  the  edges  are  very  far  away  from  the  electrodes  and  active  region;  but  serves  to 
keep  all  charged  particles  inside  the  system.  No  dielectrics  or  electrodes  have  been 
introduced  in  this  set-up.  A  brief  search  of  numerical  matrix  techniques  did  not  reveal  a 
unique  direct  (implicit)  method  for  solving  this  pentadiagonal  matrix.  It  could  however 
be  solved  using  a  computationally  intensive  Gaussian  elimination  scheme.  It  would  be 
interesting  to  actually  compare  the  time  for  convergence  of  the  iterative  line  method  used 
in  this  project  and  the  time  for  a  solution  to  the  pentadiagonal  system  using  a  direct 
method,  or  other  type  of  iterative  technique. 
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Figure  40.  Pentadiagonal  matrix  for  implicit  potential  solution. 
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Given  the  promise  of  DBD  device  technology  to  radically  alter  the  aviation  world 
of  today  and  the  specific  benefits  associated  with  its  development  to  the  United  States  Air 
Force,  should  compels  the  aerospace  industry  and  academic  community  to  seriously  look 
at  the  technology  as  a  way  of  preserving  America’s  air  supremacy  in  the  21st  century. 
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